/*   Adapted from Jonathan Katz's: uk.inc.model1p.prg,v 1.12 1997/08/07 
     Jason Wittenberg. 12/20/2001, witty@polisci.wisc.edu  
*/

/*** Set up Gauss ***/
new;
library cml,lpack,comp2,probs;

#include gauss.ext; 
#include cml.ext;
#include comp2.ext;

cmlset;

rndseed 123456789;            /* Fix seed */
setprior=1;                   /* 1=Use KK's prior, otherwise 0 */
_IncRun_M = 500;              /* This is MSIMS- set to 500 */

"Number of simulations: " _IncRun_M;
" Using KK's priors: " setprior;


/* Format stuff */
fmtt;

/*** Set up the data ****/

#include loaduk.prg;

origpca=pca;                   /* Save district ID for original dataset */

/** Intialize some matricies for storage **/
declare matrix ttwev1p,ttwey1p,ttwev2p,ttwey2p,ttwev3p,pyr,pcap,pcpp;
declare matrix cov5992p,ConEff5992p,LabEff5992p,LibEff5992p;
theta5992p = {};
ttwev1p = {};
ttwpv1p = {};
ttwev2p = {};
ttwpv2p = {};
ttwev3p = {};
ttwpv3p = {};

ttwey1p = {};
ttwey2p = {};

pyr = {};
pcap = {};

pcpp = {};
msep = {};

x1p = {};
v1p = {};
x2p = {};
v2p = {};
v3p = {};
  
tdfp = {};
ts11p = {};
ts12p = {};
ts22p = {};

fullyp = {};

conmean = {};
labmean = {};
libmean = {};
constd = {};
labstd = {};
libstd = {};

Overmean = {};
Overstd = {};
Overci = {};

posteff = {};                    /* Will hold posteriors of Coneff, Labeff, Libeff */


dfpost = {};                     /* This will hold posteriors */

P=2;
NC=0.001;

sammarkp=0;

"Start time: " time;

/** Log File **/

output file = model.out reset;

/*** Loop over elections ***/
for i (2,11,1);                  /* Should be 2,11,1 for all elecs */

  /** Set up data for this run **/
  con1  = varget(vrs[i-1,1]);
  lab1  = varget(vrs[i-1,2]);
  lib1  = varget(vrs[i-1,3]);
  con2  = varget(vrs[i,1]);
  lab2  = varget(vrs[i,2]);
  lib2  = varget(vrs[i,3]);
  ninc2 = varget(vrs[i,5]);
  pca   = origpca;

  /** Set Estimation Options **/
  _IncRun_unc = .15;
  _IncRun_NC = 1e-3;

  _IncRun_Verbose = 0;
  _RegMT_truncate = 1;
  _RegMT_psig1    = { -1.7223, .20    };
  _RegMT_psig2    = { -1.5703, .24 };
  _RegMT_pdf      = {  1.9665, .60 };
  _RegMT_prho     = {  1.3601, .35 };

  _RegMT_pb1      = {  0.0000 1e+5,  /* Constant */
                       0.0000 1e+5,  /* YL1      */
               0.0000 1e+5,  /* YL2      */
               0.0400 .008,  /* ConInc   */
               0.0100 .008,  /* LabInc   */
              -0.1500 .080 }; /* LibInc   */

  _RegMT_pb2      = {  0.0000 1e+5,  /* Constant */
                       0.0000 1e+5,  /* YL1      */
               0.0000 1e+5,  /* YL2      */
               0.0200 .008,  /* ConInc   */
               0.0600 .008,  /* LabInc   */
              -0.2500 .080 }; /* LibInc   */

 
  _RegMT_title = "Results for "$+vrs[i,6];

  _cml_CovPar = 1;

  /** Run Estimation **/
  
  { theta,cov,coneff,labeff,libeff,ttwpv,ttwev,yhat,pcp,mse,dta } =
              incrun(con1,lab1,lib1,con2,lab2,lib2,ninc2,pca);  
              
  
  nyr=rows(ttwev);
   
  /* Collect usefule quantities from each run */
            
  ttwev1p = ttwev1p | ttwev[.,1];
  ttwpv1p = ttwpv1p | ttwpv[.,1];
  ttwev2p = ttwev2p | ttwev[.,2];
  ttwpv2p = ttwpv2p | ttwpv[.,2];
  ttwev3p = ttwev3p | ttwev[.,3];
  ttwpv3p = ttwpv3p | ttwpv[.,3];
  
  ttwey1p = ttwey1p | yhat[.,1]; 
  ttwey2p = ttwey2p | yhat[.,2];
 
  pyr   = pyr | years[i].*ones(nyr,1);
  
  pca = lget(dta,3);
  
  
  pcap = pcap | pca;
  
  pcpp = pcpp | pcp.*ones(nyr,1);
  msep = msep | mse.*ones(nyr,1);
  
  
 tdfp = tdfp | tdfi(P,theta).*ones(nyr,1); 
 Sigma = tsigi(p,theta);
 ts11=sigma[1,1];
 ts22=sigma[2,2];
 ts12=sigma[1,2];
 
 ts11p = ts11p | ts11.*ones(nyr,1);
 ts22p = ts22p | ts22.*ones(nyr,1);
 ts12p = ts12p | ts12.*ones(nyr,1);
 

  
  temp = lget(dta,1);
  x1 = temp[.,2:cols(temp)];
  v1 = temp[.,1];
  
  x1p = x1p | x1;
  v1p = v1p | v1;
  
  temp = lget(dta,2);
  x2 = temp[.,2:cols(temp)];
  v2 = temp[.,1];
  
  x2p = x2p | x2;
  v2p = v2p | v2;

  v3 = 1-v1-v2;
  v3p = v3p | v3;
 
  theta5992p = theta5992p | theta;
  
  MP = -9*ones(rows(v1),1);
  MP = substute(MP,v1.>NC .and v2.>NC .and v3.>NC,0);
  MP = substute(MP,v1.<=NC .and v2.>NC .and v3.>NC,1);
  MP = substute(MP,v1.>NC .and v2.<=NC .and v3.>NC,2);
  MP = substute(MP,v1.>NC .and v2.>NC .and v3.<=NC,3);
  MP = substute(MP,v1.>NC .and v2.<=NC .and v3.<=NC,4);
  MP = substute(MP,v1.<=NC .and v2.>NC .and v3.<=NC,5);
  MP = substute(MP,v1.<=NC .and v2.<=NC .and v3.>NC,6);

fully=(MP.==0);
fully=substute(fully,(MP./=0 .and MP./=-9),0);
 
fullyp = fullyp | fully;
 
  
  /** Store Results **/
  theta5992p = theta5992p | theta;
  cov5992p = lput(cov5992p,cov,i-1);
  ConEff5992p = lput(ConEff5992p,ConEff,i-1);
  LabEff5992p = lput(LabEff5992p,LabEff,i-1);
  LibEff5992p = lput(LibEff5992p,LibEff,i-1);
   
 /* This is taken from results1p.prg */
  
  /* Conservative Effect Year by Year */
  ConMean = ConMean | meanc(packr(ConEff))'.*ones(nyr,1);
  ConStd = ConStd | stdc(packr(ConEff))'.*ones(nyr,1);
  
  /* Labour Effect Year by Year */
  LabMean = LabMean | meanc(packr(LabEff))'.*ones(nyr,1);
  LabStd = labStd | stdc(packr(labEff))'.*ones(nyr,1);
  
  
  /* Liberal Effect Year by Year */
  LibMean = LibMean | meanc(packr(LibEff))'.*ones(nyr,1);
  LibStd = LibStd | stdc(packr(LibEff))'.*ones(nyr,1);


  /* Pooled Effect Year by Year */
  IncEff = ConEff[.,1]~LabEff[.,2]~LibEff[.,3];
  IncEff = sortc(meanc(packr(IncEff)'),1);
  OverMean = OverMean | meanc(IncEff)'.*ones(nyr,1);
  OverStd = OverStd | stdc(IncEff)'.*ones(nyr,1);
  M = rows(IncEff);
  OverCI = OverCI | (IncEff[round(0.05*M)]~IncEff[round(0.95*M)]).*ones(nyr,1);  
   
  posteff = posteff | (years[i].*ones(rows(coneff),1))~(coneff~labeff~libeff);
  
  "End time for election: " i ": " time;

endfor;

y1_1p=x1p[.,2];
y1_2p=x1p[.,3];
ConIncP=x1p[.,4];
LabIncP=x1p[.,5];
LibIncP=x1p[.,6];

/* The following variables will be saved and exported */

ds=pcap~pyr~y1_1p~y1_2p~
   ConIncP~LabIncP~LibIncP~v1p~v2p~v3p~
   ttwev1p~ttwpv1p~ttwev2p~ttwpv2p~ttwev3p~ttwpv3p~tdfp~
   ttwey1p~ttwey2p~ts11p~ts12p~ts22p~msep~pcpp~fullyp;

/** Close Log **/
 output off;
 
 screen off;
 output file=ukout1.out reset;
   ds;
 output off;
 
output file=posteff.out reset;
   posteff;
output off;
screen on;

"The program has finished. Data written to ukout1.out";
"Run IMPUTEY.PRG";





 
